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Metal - Insulator Transition in 3D Quantum Percolation 



We present the metal - insulator transition study of a quantum site percolation model on 
simple cubic lattice. Transfer matrix method is used to calculate transport properties - 
Landauer conductance - for the binary distribution of energies. We calculate the mobility 
edge in disorder (ratio of insulating sites) - energy plane in detail and we find the extremal 
critical disorder somewhat closer to the classical percolation threshold, than formerly 
reported. We calculate the critical exponent v along the mobility edge and find it constant 
and equal to the one of 3D Anderson model, confirming common universality class. 
Possible exception is the center of the conduction band, where either the single parameter 
scaling is not valid anymore, or finite size effects are immense. One of the reasons for 
such statement is the difference between results from arithmetic and geometric averaging 
of conductance at special energies. Only the geometric mean yields zero critical disorder 
in band center, which was theoretically predicted. 

Keywords: Disorder, mobility edge, critical exponent. 

1. Introduction 

Transport properties of disor dered solids attracted lot of attention. The quantum 
percolation (QP) problerrP^ has been studied several dec ades and some questions 
are still open, compared to the more settled Anderson modePEl. First of all, do both 
models belong to the same universality class? Positive answer means that the critical 
exponent v of localization length should be the same, i. e. va — 1.57(2) in 3D and in 
presence of time reversal symmetry. The best data come from transfer matrix (TM) 



parameter scaling (SPS). It is possible, that SPS does not work in the mid dle o f the 
conduction band for 3D QP, which is the case also for ID Anderson modef^, but 
not for 3D Anderson model (at least numerically). Contrary to the still discussed 
2D QP case^, all authors agree that metal - insulator transition (MIT) is present 
in 3D, but they often disagree in details. 
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The mobility edge trajectory (MIT li ne in energ y E vs. disorder plane) in 3D 
is qualitatively clear and widely accepted I 9 | 10 | H I 121 w ith possible jumps at special 
energies of incoming non-interacting transport careers (E = 0, ±1, ±\/2, ...)r^H We 
aim to comment the nature of these jumps in more detail - for site percolation, but 
the bond percolation behaves similarly!^. 

The common universality further requires constant v along the mobility edge. 
The exponent shoul d be the same both with respect to changing E at fixed disorder 
and vice verscP3HSl. We will check all these statements numerically. Another widely 
accepted fact is that v should be identical both for site and bond percolation. 

The values of the exponent v given in literature differ heavily. Some methods 
seem to overestimate the true value, if we accept va as correct also for QP. Renorm- 
group calculations gave v = 2.2 ^ or v = 1.81(6) Thouless conductivit}0^ 
yielded v = 1.95(12) on fee lattice. On the other hand a special version of TM 
method with approximative reduction of matriceJ^ suggested v = 0.75(10) and 
cluster expansion with Pade approximation^^ yielded v = 0.38(7). The last value 
violates even the Chayes bound v > 2 /d, where d = 3 is dimension. It is known from 
the Anderson model, that level spacing methods tend t o so mewhat underes tim ate 
v compared to the TM value; they suggest v = 1.35(10) ^] and v = 1.45(7) ^ for 
3D QP. Thus we can say that the latter results confirm the common universality. 
The same is true of calculations using retarded Green's function method^, with 
resulting v — 1.6(2) at E — 2, but for small system sizes L < 8 and rather spread 
data. Another paper of the same groupM^ reported lower v — 1.2(2). The main goal 
of this paper is to present another, TM based confirmation of the assumption, that 
Anderson model and QP belong to the same universality class. 



2. Model and method 

Let us recall the tight-binding Hamiltonian: 

H = Y / e i \i><i\+ Uj\iXj\ (1) 

i <i,j> 

where the first sum counts the L 3 lattice sites, the second sum goes over near- 
est neighbors and we set the hopping t^ = 1, thus fixing the energy scale. The 
length scale is set by choosing the lattice constant a = 1. For a binary system the 
distribution of on-site energies €i has the form: 

P{t t ) = (1 - q)8{e l - e A ) + qd(e z - eg). (2) 

We prefer to use q, the ratio of insulating atoms B, instead of often used p = 1 — q 
because q is the measure of disorder, induced among perfectly conducting A atoms. 
The quantum percolation model can be got either by removing appropriate ty- — * 
in Eq. (TTJ) or by sending es — * oo. The latter case is more simple for TM method, as 
we need all ty = 1 in the transfer direction and it suits better for site percolation. 
Thus we set the on-site energy for A-type atoms to €a = and the one for the 
randomly positioned B-type atoms to a large number^H], typically es = 10 7 . The 
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prize one has to pay in numerical approach are ill-defined matrices and a non-zero, 
but very small tunelling conductance through B sites with still finite barrier - though 
the latter disadvantage shows up only in deep insulators^, behind the scope of this 
paper. 

For Anderson model the P(ei) in Eq. ((2|) is continuous, e. g. box distribution 
from an interval of the width W (disorder): e, G [-W/2, W/2]. It should be further 
mentioned, that both models are fully symmetric with respect to E — > — E, thus we 
will restrict our investigation to one (left) half of the band: —6 < E < 0. 

The classical site percolation threshold for simple cubic lattice is qc — 0.6884 
and it is known, that the quantum critical value qq < qc for the whole mobility 
edge. The quantum wave requires a better highway than the narrow p ath of the 
percolation cluster backbone. Up to now the maximum calculated 

value™ 

was 

qq. 0.56(1), but we will show that the difference from qc is smaller. 

Our numerical calculations are based on the usual transfer matrix method, de- 
scribed elsewhere^. The conductance g in units e 2 /h is given by the Landauer 
formula g = 2 Tr t + t. Here t is the 1? x L 2 transmission matrix, specifying the 
transmitted (not reflected) part of L 2 planar waves, entering the sample on the 
left and leaving it on the right, t is extracted from a 2L 2 x 2L 2 transfer matrix. 
It was mentioned above, that we have to do with ill-defined matrices, even more 
than within Anderson model. In transfer matrix approach we deal rather with t -1 
matrices and the precision of their small eigenvalues gets quickly lost - but we need 
them accurate, as they become important after recovering the t matrix. An efficient 
way how to treat this problem was described in Refs. [24] and [4j we have to perform 
the renormalization procedure at least after each two TM multiplications, but then 
we are able to include larger samples up to I = 20. This is useful, as finite size 
effects turn out larger than for Anderson model. The most numerically sensitive 
region is around the band center, i. e. for small \E\, where we have to perform the 
procedure after each single TM multiplication^. 

Once we have the matrix tH, we can calculate also its eigenvalues r^, satisfying 

N N ^ 

where the positive quantities Xi play similar role for cubic samples, as Lyapunov 
exponents do in quasi-lD geometry. The smallest one is labeled X\. 

In localized samples, the Landauer conductance decays quickly with sample size, 
g ~ exp(— L/£). The correlation length £(q,E) diverges at critical point (Eq,pq) 
as £ ~ (q — qq) ~ v for Eq fixed and vice versa. The SPS assumption implies one 
combined variable from (q, L) or (E, L) and in the vicinity of the MIT point we 
expect g(q,L) ~ g(L/£(q)) ~ g(L (q — qo) v ). If we perform Taylor expansion in 
(<7 ~ Qq) variable and add some finite-size correction, we get 



g = g Q + A x {q - q Q )L^ + A y L~ v + o([(g - qg)^} 2 ) , (4) 
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Fig. 1. Conductance g(es) near special point E = —1. Two samples: the first, delocalized one 
- overlaying full symbols above; another one - empty symbols reaching from slightly to deeply 
localized case. 

where y is an irrelevant exponent and its term becomes negligible for large L. 
Analogous relations can be expected for q <-> E, i. e. with the single variable (E — 
Eq)L x / v and qQ fixed. The conductance g will be mostly an ensemble average 
(g), taken over many samples, typically 20.000. It is knowrP^ that many other 
quantities would fit Eq. (2]), e. g. the geometric mean exp((ln g)),l/(l/g), percentiles 
of the distribution of conductances P(g), (xi), current momenta^ r"), inverse 
participation ratioi^, etcff^l Nevertheless, quantities like (In g) or (xi) are more 
sensitive to numerical instabilities, especially for larger p and L or very small \E\. 

3. Numerical results 

It was showrff^l, that the density of states (DOS) has sharp peaks at special ener- 
gies, \E\ = 0,1, \/2, (3 ± v / 5)/2, ... These eigenvalues belong to submatrices of the 
Hamiltonian Eq. (fT]), created by small islands of A sites completely surrounded by 
B sites, i. e. isolated from the rest. It was found, that g can quickly fall down at 
these energies and they can influence the shape of mobility edge. 

Before going to main results, let us check our numerics in the vicinity of the 
special point E = — 1. Conductance g should be cb independent for large enough 
cb- We can see in Fig. [TJ that for delocalized sample with g > 1, common upper 
line, this is true for €b > 10 4 even at E = —1. For a localized sample, all lower lines, 
the value depends strongly on the energy E close to the special E = — 1, but still a 
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plateau develops for larger es- Exception is the special point itself, where g oc , 
in other words g — ► for es — * oo. The same g(£B) dependence is found for a 
test sample with one perpendicular square barrier of B sites completely stopping 
the transport through otherwise perfect cube of A sites; thus this dependence is 
an unusual form of tunelling. The picture is qualitatively the same at E = 0, just 
for much lower disorder q. It is worth mention, that in the latter case we sur ely 
do not have to do with simple quantum tunelling through closed finite samples^ 1 
(closed in the sense of classical percolation, regardless of E), as we chose samples 
with q — 0.37, far bellow qc- Alternatively formulated, the same samples become 
metallic for larger \E\. One can also see that the geometric mean of conductance 
tends to zero though the arithmetic (<?) remains positive at special points. We chose 
also E = —y/2 and after some effort we found a localized sample of the type as in 
Fig. [T] bottom lines, but the ratio of such specially localized samples becomes much 
lower than for E = 0, ±1. 

The possible deep localization at special energies was explained in several ways. 
In Ref. [T] it was supposed, that the clusters of insulating atoms around con- 
ducting islands could have extremely large reflectivity at these energies. Another 
explanation^ argued that localization on dead ends of the spanning cluster is strong 
at special energies. Let us add yet another possibility. The dispersion law gives 

ima irma . , . 

E = — 2t cos k z a — 2t cos 2t cos , < k n < n/a; n, m = 1, .... L (5) 

L+l L+l -«-/'' . > w 

where k z is the wavenumber in transfer direction z and the wavenumbers in perpen- 
dicular directions x, y are quantized by hardwall boundary conditions. For special 
values of energies the whole vector (k x , k y ,k z ) can point in (body) diagonal direc- 
tion and the planar wave is then fully reflected by insulating site, siting exactly 
in this direc tion. If all planar waves are fully reflected, the sample becomes deeply 
localized^. 

Summarizing, with cb — 10 7 , we always get correct g for delocalized samples and 
also for medium localization. Strongly localized samples require special attention. 
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Fig. 3. The mobility edge, or MIT line. 

Now let us find one critical point in detail. For q < qq most samples are metallic 
and (g) oc L d ~ 2 oc L, for q > qQ most samples are localized and (g) oc exp(— £/£)• 
At MIT, g becomes L independent. We chose the value qQ = 0.57, and found 
Eq = —2.36 from the crossing point in Fig. and [2b)- As qQ is far enough from 
qc, and Eq is not close to any important special value (in the sense of relevant 
peaks in DOS), we can use even the (xi) as order parameter. The data were fitted 
according to Eq. fH and the resulting exponent is v = 1.6(1) for the (q — qQ) 
dependence and v = 1.5(1) for the (E — Eq) one, both in perfect accordance with 
va- Of course, taking (g) instead of (xx) yields the same critical point and v. 



Having performed these calculations for almost two dozens of fixed energies, 
typically changing only q, we can plot the mobility edge in Fig. [3l Some pairs of 
numerical values can be found in Fig. [6l The point (—6, 0) was added artificially - 
there are no open channels there. It was already mentioned above that the situation 
is mirror symmetric in E — » — E, we do not plot the obvious interval < E < 6. 
The upper border of the graph was set to the classical qc- The maximum v alue of 
our MIT line is cca q — 0.603(2), greater than 0.56(1) value reported before™, 
but with smaller samples L < 9 and L < 8, resp. 

Now we will compare our results to those of Ref. [T3] They introduced the local 
DOS for each lattice site and then performed two averages at once (in two manners). 
The first one was ensemble average as we do, the second one was averaging with 
respect to all sites. The two manners were geometric and arithmetic mean. At last 
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Fig. 4. (a) Detailed view of the mobility edge at band center, (b) L - independent crossing point 
of arithmetic mean (g(p)) at band center. 



they calculate the ratio of geometric/ arithmetic mean and plot lines, where this 
order parameter is constant. The line entering the point (—6, 0) compares very 
well with our mobility edge in Fig. [3l with their minimum at cca p — 0.4, i. e. 
our maximum at q = 0.6 and a rather wide plateau around this maximum. The 
main difference is a discrete set of energies, where their order parameter and thus 
also the critical line drops to zero. But the only cases, where they could show 
this explicitly where integers E = 0, ±1. With irrational values like E = ±V2 the 
localized case is hard to find numerically, as its region is extremely narrow in energies 
and also the corresponding peaks in DOS are much smalleJISl than those at three 
integer E. Hence we plot the detailed view of the band center in Fig. a), where 
both calculations of critical qQ with the use of geometric and arithmetic mean of 
conductance g are presented - of course only ensemble averaged. The critical disorder 
qQ from calculations with geometric mean of g is systematically somewhat larger 
outside the band center, but only by several percent - beware that the q axis is 
not drawn to zero. For very small \E\ = 10~ 6 , 10~ 7 the qQ from geometric mean 
already becomes slightly smaller than that from arithmetic mean and directly at 
E = it drops steeply to zero. It should be stressed out, that the decay of MIT line 
at band center down to cca qQ = 0.37 is independent of the type of averaging, it 
happens within a clearly non-zero interval of energies and it should be distinguished 
from extremely narrow drops of MIT line to zero at special energies, present only 
with geometric averaging. Similar relatively broad interval of energies, where the 
mobility edge should decay, was reported^ around E — 1 . With arithmetic mean 
and €b = 10 7 , we do not see it at all. We should also say that the standard mean 
DOS, calculated up to L = 16, has clearly non-zero value for q = 0.37 in the 
vicinity of the band center, apart from high central peak at E = 0. According to 
our calculations, only at cca q — 0.44 the DOS drops to zero for small \E\ and a 
gap opens for even larger q. 

The exponent v remains constant along the whole mobility edge line, just the 
error bars rise as we approach the borders of the band \E\ = 6 and E = 0. But even 
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Fig. 5. Critical values of (g), (lng) and (si) along the mobility edge. 

for E = —0.001 and = 0.39 it was still v = 1.5(3). Important exception is the 
point E = itself. Only with geometric mean, it has q<Q = 0, as suggested already 
in Ref. [TJ One possible explanation is that single parameter scaling (SPS) does not 
work here anymore. The argument is, that geometric mean of g differs substantially 
from the arithmetic one. Furthermore, if we insists on SPS and arithmetic (g) as 
order parameter, standard analysis of data from Fig.^b) according to Eq. (J4| yields 
qq = 0.37 at E = and v w 1 with (q — <7q) variable. We tried also to calculate v 
with (E— 0) variable, but the g(E) dependences changed their slope below \E\ < 10 7 
and the analysis was too unstable. It is improbable, that common universality would 
be broken at one spectral point, the v from Fig. 0|b) data smaller than va could 
rather mean, that SPS is broken at E — 0. But one can also argue, that finite size 
effects are so strong here, that even sample sizes up to L — 20 are insufficient. This 
would be exceptional, though not completely excluded. In Ref. |T5] the authors had 
similar problems with 3D Anderson model at the band edge (close to E = 6) and 
explained them by insufficient available sample sizes. Summarizing, we always get 
a well defined crossing point of mobility edge, where data become L-independent 
as they should at criticality, even for E = 0, see Fig. Q|b). Just the slopes of the 
crossing lines behave non-universally in band center, at least for available system 
sizes. 



Let us present mean critical values of several quantities, namely (g)c, ( m .9)c 
and {x\)c- For g < 2 we can roughly say that it is given mainly by the largest 
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eigenvalue g ~ 1t\ = 2cosh _2 (xi/2) and after simple algebra lng w 3 In 2 — x±, 
which is well fulfilled in Fig. [5j of course only outside the band center. For small 
\E\ we get gc > 2 as a sum of several important T\. 

It should be emphasized, that (g)c remains finite, but (lng)c and {x\)c diverge, 
as we approach E — 0, see comments to Fig.[TJ The divergence is given by a minority 
of fully localized samples with g — > for cb — > oo , which do not really influence the 
arithmetically averaged gQ ax ~ 6.5. 

Further we present several critical conductance distributions Pc(g) and PcO^g) 
along the mobility edge. The distribution at E = and p = 0.37 was added for 
comparison, too. 



0.8 




Fig. 6. Critical distributions of g (a) and lng (b) along the MIT line. 



We do not plot critical distributions for —0.05 > E > —3 which are mutually 
almost overlaying, because of similar qQ . Further we can see typical non-analyticity^ 
at g = 2 or lng = In 2, if P(g) or P(lng), resp., is already decaying there. The 
Pc(lng) at E = —3 and those with maximal qQ look similar to Pc(\ng) of 3D 
Anderson model, but they are not identical for any (E, q) - the SPS works separately 
inside the models, not universally. This is known also for anisotropic Ande rson 
models and those in other dimension than d — 3. Even for the bifractaP^SI w ith 
the lowest d = 2.226 and gc ~ 7 the Pc{g) looked more Gaussian than those in 
Fig(6|a) on the right. This is connected with far larger var g = (g 2 ) — (g) 2 of 3D 
QP compared to low - dimensional Anderson model at criticality or to universal 
conductance fluctuations in metallic regime. 

Let us shortly discuss the metallic regime. Away from band center, say with 
E = —1.05, we found typical metallic behavior for medium disorder 0.2 < q < 0.5, 
i. e. P(g) practically Gaussian, (g) oc L and var g constant (L and q independent), 
with values cca 10 percent above the theoretically predicted universal conductance 
fluctuations^. Even smaller q brings traces of ballistic behavior, var g becomes 
system-size dependent (rising with L). Contrary to these standard properties, at 
the band center £ = 0we did not find typical metal for any q and available L. We 
still have (g) oc L, but now var g oc L instead of reaching the universal constant 
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and cca one order of magnitude greater - rather like the mentioned almost-ballistic 
behavior (perfect ballistics would have (g) oc L 2 ). Also P(g) is somewhat deformed, 
not Gaussian. The transition from this behavior to standard metal at larger is 
continuous, var g decays with rising \E\. If this is only finite size effect, than it is 
again immense at band center. 

4. Conclusions 

We already stated recently 21 , that 2D quantum percolation problem and 2D An- 
derson model belong to the same universality class, meaning that d — 2 is lower 
critical dimension for both models and there is no MIT in time reversal symmetry 
case. Here we add the confirmation of common universality in 3D, implying identi- 
cal critical exponent which we have shown within numerical accuracy along the 
mobility edge, with possible exception of the conduction band center E = 0. An- 
other new result confirms equal v when cutting the mobility edge (MIT line) either 
in disorder or in energy direction. All these properties result from single parameter 
scaling hypothesis and they get lost at band center, where either the SPS does not 
work anymore, or we face immense finite-size effects, that cannot be overcome with 
system sizes up to L — 20. This statement is supported also by loosing many typ- 
ical properties of metallic samples at band center. Concerning the mobility edge, 
we found the maximum critical disorder qq slightly above 0.6, which is more than 
q Q w 0.56 reported 

before™]. 

The shape of the mobility edge at band center de- 
pends on the choice of order parameter. The arithmetically averaged (g) (commonly 
preferred) gives clearly non-zero qq f=a 0.37 at E = , whereas the geometrical mean 
gives qq = in accordance with other numericA^ and theory^. We present a new 
alternative explanation of extremely low conductance for a minority of samples, 
which push the geometric mean to zero. It is based on the fact, that for special en- 
ergies the wavevectors of incoming planar waves point in body diagonal directions 
and the waves can face exactly the insulating sites, followed by full reflection. 

The rational special energies E = 0,±1, have clearly different arithmetic and 
geometric means and also the resulting qq. The irrational ones, e. g. E = ±72, 

seem 

numerically almost insensitive of the type of averaging, as the samples with nearly 
zero conductance become very rare. But the rather steep descent of mobility edge 
at band center takes place at non-zero interval of energies with SPS still working for 
both types of averaging and it should be distinguished from the extremely narrow 
regions at all special energies, apparent only for geometric mean of conductances. 
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